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ABSTRACT 

The made-to-measure A^-body method (Syer & Tremaine 1996) slowly adapts the particle 
weights of an A^-body model, whilst integrating the trajectories in an assumed static potential, 
until some constraints are satisfied, such as optimal fits to observational data. 1 propose a novel 
technique for this adaption procedure, which overcomes several limitations and shortcomings 
of the original method. The capability of the new technique is demonstrated by generating 
realistic A^-body equilibrium models for dark-matter haloes with prescribed density profile, 
triaxial shape, and slowly outwardly growing radial velocity anisoti'opy. 
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1 INTRODUCTION 

A standard problem in contemporary galaxy dynamics is the in- 
terpretation of kinematic observations of galaxies in terms of their 
orbital structure as well as their dark and luminous matter distri- 
bution. There are several methods one can employ for this prob- 
lem. First, moment-based methods find solutions of the Jeans equa- 
tions (or higher-order velocity moments of the coUisionless Boltz- 
mann equation) that best fit the observed moments, such as den- 
sity and velocity dispersion. Secondly, distribution-function-based 
methods directly fit the distribution function to the data, which 
can be more general than mere moments, e.g. the line-of-sight ve- 
locities of many individual objects. Both techniques are usually 
restricted to spherical or, under certain simplifying assumptions, 
axisymmetric systems (though for different reasons^). However, 
distribution-function models are technically much more challeng- 
ing (since an integral equation has to be solved instead of differen- 
tial equations) and hence much less used than the moment-based 
approach. In both cases, astrophysically unjustified assumptions, 
such as velocity isotropy, are often made in order to make the prob- 
lem tractable. 

Thirdly, Schwarzschild's (1979, 1993) orbit-based method 
constructs a dynamical model by first integrating many orbits over 
many orbital times in an assumed gravitational potential, whereby 
recording their properties in an orbit library, and then superposing 
them such that a best fit to the data is obtained. This is a pow- 



* Email: walter.dehnen@astro.le.ac.uk 

^ For moment-based models, symmetry reduces the number of indepen- 
dent moments and enables simple assumptions necessary to close the Jeans 
equations. For non-spherical distribution-function models, knowledge of 
isolating integrals of motion other than energy is necessary, and angular 
momentum is available only for axial symmetry. 



erful method, since it comes, in principle, without restrictions on 
the symmetry, and one may even obtain the distribution function 
(Hafner et al. 2000). However, in practice most applications are 
restricted to axisymmetry, since there are several technical sub- 
tleties to overcome when applying the method to potentials with 
a complex phase-space structure, the typical situation for triaxial or 
barred systems (though this is not impossible and has been done, 
e.g. Hafner et al. 2000). 

Fourthly, in 1996, Syer & Tremaine (hereafter ST96) in- 
troduced the 'made-to-measure A'^-body method', which slowly 
adapts a TV-body model to fit the data. The first application of this 
method came as late as 2004, when Bissantz, Debattista & Ger- 
hard used it to construct a dynamical model for the Milky Way's 
barred bulge and inner disk. More recently, De Lorenzi et al. (2007, 
hereafter DL07) refined the method to incorporate observational 
errors; this has since been applied for modelling elliptical galax- 
ies to asses their dark-matter content (De Lorenzi et al. 2008a,b). 
The made-to-measure (hereafter 'M2M') method is as powerful as 
Schwarzschild's orbit-based method, and in fact is closely related. 
Whereas in Schwarzschild's approach orbits are first separately in- 
tegrated and then superimposed, these two steps are merged in the 
M2M method: trajectories are integrated and their weights adapted 
at the same time. As a consequence there is no need for an orbit li- 
brary and all the technical difficulties associated with it. However, 
with the M2M method as proposed by ST96 & DL07 some prob- 
lems remain, as I shall discuss, in particular the appropriate time 
scale for adapting the particle weights of the A'^-body model. 

Finally, Rodionov, Athanassoula & Sotnikova (2009) intro- 
duced a variation of the M2M technique (though the authors did not 
make this association), which they dubbed the 'iterative method'. 
Their method starts from a near-equilibrium dynamical model (con- 
structed by any other method), which is alternately relaxed under 
self-gravity (to evolve towards equilibrium) and adapted to pre- 
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scribed properties. In practice, tliis metiiod too employs the A'^- 
body approach and, like the traditional M2M technique, suffers 
from the time-scale problem. 

Another application of all the aforementioned techniques is 
the generation of A'^-body initial conditions representing a galaxy 
or galaxy component (though moment-based models additionally 
require the incorrect assiunption of a Gaussian velocity distribution 
and should not be used for this purpose). However as mentioned 
above, distribution-function models, which are the most commonly 
employed technique for generating A^'-body galaxy models, are re- 
stricted to spherical (or, under simplifying assumptions, axial) sym- 
metry, which seriously limits their realism. Since the M2M tech- 
nique works directly with A'^-body data, it offers a simple and natu- 
ral way to generate A^-body initial conditions with prescribed prop- 
erties (ST96), in particular non-spherical shape and non-isotropic 
velocity structure. As we shall see, however, this requires some 
modification to the traditional M2M technique. 

In this paper, I revisit the M2M method and propose several 
modifications aimed at improving it, in particular in view of its ap- 
plication for tailoring A^-body initial conditions. I present the tradi- 
tional ST96 & DL07 version of M2M in section 2 and my modifica- 
tions to the method in §3, while §4 presents some tests of tailoring 
non-spherical and/or velocity-anisotropic AT-body models. Finally 
§5 discusses the results and concludes. 



2 TRADITIONAL M2M 

In this section, the M2M method as laid out by DL07 (which in 
turn was based on ST96) is outlined, though with slightly different 
notation and conventions. 

The fitting of the A^-body model to the data is expressed 
as maximisation problem: the AT-body model shall maximise the 
merit function 



he. 



(1) 



Here, C is the constraint function, which measures the goodness of 
fit of the AT-body model to some target. There are many possible 
choices for C, but for now let us follow DL07 and consider a x^- 
hke measure of the deviation of moments of the AT-body model 
from target values 



Here, 



(2) 



(3) 



are some moments of the model defined via the kernel Kj{xi,Vi) 
and the particle weights Wi = rrii/Mtot, while Yj are the targets 
for those moments and represent the observed data^ with uncertain- 
ties aj . 

Simply minimising C is not a well defined procedure for two 
reasons: first, there is no point in reducing C well below the ex- 
pectation value even if this were possible (this would amount to 

^ This standard practice restricts the data to be just moments of the dis- 
tribution function, and excludes, for example, the line-of-sight mean and 
dispersion velocity, which are functions of moments. However, this restric- 
tion is not fundamental and the method can easily be extended to fit any 
function of moments, see §4.2 for an example 



'fitting the noise'); second, minimising C may not be uniquely con- 
straining the AT-body model: there are, for instance, many possible 
equiUbrium models with the same density. Thus, in order to yield a 
well-defined problem, one has to regularise the merit function by a 
penalty functional S times a Lagrange multiplier /i, which controls 
the amount of regularisation. The penalty function is traditionally 
taken to be the pseudo-entropy 



5 = -^<log^ 



(4) 



with Wi = Wi/J2j '^3 the normalised weights and {wi} a pre- 
determined set of normalised weights, the so-called priors. For 
general priors, S defined in this way is the Kullback & Leibler 
(1951) information distance (also known as 'K-L divergence') of 
the model corresponding to Wi = m from the actual AT-body 
model, i.e. S penalises against deviations of the normalised weights 
from the priors. Only for wt oc f^'', where fi denotes the value of 
the equilibrium distribution function corresponding to Wi = 1 /N, 
does S reduce to the true entropy of the AT-body model (plus a 
constant; this corrects statements made by ST96 & DL07). 

The idea of the M2M method is now to adjust the weights 
slowly such that Q is maximised. The standard method is to evolve 
the weights according to 



Wi = ewi Ui 

with some rate of change e and the velocity of change^ 
dQ 



(5) 



U, 



dwi 



(6) 



For the particular choice (2) of the constraint function, this gives 

dS ^Y^-_yj 

.2 

j=l -3 



dS Yi — y-j , 



(7) 



For sufficiently small e, integrating (5) will increase Q and even- 
tually result in a AT-body model for which Q is maximal and the 
Wi no longer change. This method is similar to (and was in fact in- 
spired by) Richardson (1972) - Lucy (1974) iteration, though with 
a much reduced step size. 

Unfortunately, it is not as simple as that, because the merit 
function, being a function of the randomly sampled particle tra- 
jectories, is itself a random variable and fluctuates even with fixed 
weights. In order to suppress these fluctuations, traditional M2M 
replaces the model moments j/j in (2) with their time-averaged val- 
ues yj, which are obtained by integrating the differential equation 



(8) 



starting with* yj = yj at t = 0. If fitting to observed data with 
finite uncertainties aj, this method has the virtue that the model 
uncertainties due to AT-body shot noise (which have been ignored 
in the definition of the constraint function) are much reduced. 

A problem with this time-averaging is that the computation of 
the derivatives required for the velocity of change (6) is no longer 



^ Unfortunately, ST96 dubbed ewiUi the 'force of change', which is an 
inaccurate analogy since it is proportional to the first time derivative of the 
dependent variable. Below I introduce a method which indeed uses the sec- 
ond time derivate, for which the expression 'force of change' is much more 
appropriate. 

Corresponding to yj (t) = yj (0) e"'*' +a /J e'^(*'-*)j/j (t') dt'; ST96 
& DL07 give y = a e"(*'~*'2/j(t') dt' which results from integrat- 
ing (8) from t = —oo, a practical impossibility. 
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straightforward. In fact, dyj/dwi is simply the time-averaged ker- 
nel — the quantity which in Schwarzschild's method is stored in 
the orbit library. In traditional M2M, one simply replaces yj with 
yj directly in equation (7). Even though, this means that weight 
adaption is not strictly along the gradient of the merit function, it 
appears that the method still converges in practice, though it is not 
completely obvious that it always does (DL07 introduce the time- 
averaging only after they argue for convergence), in particular for 
other forms of the constraint function than simple expressions 
on model moments. 



3 A NOVEL M2M METHOD 

In this section I criticise the traditional M2M method and propose 
alternatives and/or modifications, which ultimately cumulate in a 
novel method. 



While this may be a viable method, I propose a somewhat 
different approach which incorporates the total-weight constraint 
into the adjustment step. I start by observing that the unconstrained 
maximum of the modified merit fimction 

Q*{'w) = Q{w*) +\nJ2k^k -Y^f,Wk, (12) 

maximises Q{w) subject to the constraint 1 = '"'k (e g- 
Dehnen 1998). Thus, the total-weight constraint can be incorpo- 
rated by replacing Q with Q* in equation (6). Note that since Q* 
depends on the constraints only through the w*, the AT-body sys- 
tem must still be re-normaUsed after each adaption step, but the 
step hardly carries the system away from normalisation. 

Alternatively, the total weight of the A'^-body system may be 
allowed to float freely and be constrained only by the data via the 
constraint function. 



3.1 Time scales 

An important issue is the appropriate choice for the adjustment rate 
e. The velocity of change (6) varies on the orbital time scale of 
the ith trajectory, because different parts of the orbit contribute dif- 
ferently to the merit function. Clearly, the weight adaption should 
happen adiabatically, i.e. c < fij with fi, the natural (orbital) fre- 
quency of trajectory i. Unfortunately, the orbital frequencies of the 
N trajectories may easily vary by many orders of magnitude, such 
that meeting this condition for all of them becomes a serious prob- 
lem. In traditional M2M this is not really solved: using a very low 
adaption rate e ensures that the weights for all but the outermost 
trajectories are adapted adiabatically. 

One may think that using individual adjustment rates d ex. Qi 
would solve the problem. However, this is not the case: the method 
no longer converges (it does initially, but eventually convergence 
stalls well before reaching the optimum), presumably because such 
an alteration changes the direction of adjustment away from the 
gradient of Q. Instead, I turn the tables and achieve uji oc e by 
integrating each trajectory on its own dynamical time scale. To this 
end, I introduce the dimensionless time 



t/Ti 



(9) 



with Ti — 2n/Qi the orbital period, such that with respect to this 
new dimensionless time each trajectory has natural frequency uj = 
2n. The equations of motion for the A'^-body system expressed in t 
are 



(10) 



where a prime denotes derivative w.rt. t. Conversely, the M2M 
equation remains: 

w'i = eWiUi (11) 

such that now e is a dimensionless rate per orbit for each particle. 
In practice, a rough estimate for the orbital period, e.g. based on the 
epicycle approximation, is sufficient for Ti. 

3.2 Enforcing total-weight conservation 

With the traditional M2M formulation, conservation of the total 
weight is not guaranteed, as the maximum of Q may occur at 
Wi ^ 1. This problem has not been discussed by ST96 & 
DL07, and I assume that it is dealt with by simply re-normalising 
the weights after each step. 



3.3 An alternative adjustment 

As discussed in the last paragraph of §2 the traditional time- 
averaging procedure interferes with the computation of the gradi- 
ent of the merit function. Moreover, as I shall discuss in the next 
subsection, the time-averaging is particularly undesirable when us- 
ing the M2M method for tailoring AT-body initial conditions. These 
considerations lead me to consider a different time-averaging ap- 
proach: instead of averaging the moments, I consider suppressing 
fluctuations in the merit function by averaging Q* (or Q) itself and 
its derivatives. In analogy to the moment-averaging equation (8), 
this would yield 



dwi 



with the dimensionless averaging rate r). Combining this with the 
weight-adjustment equation (11) results in a second-order differen- 
tial equation for ipi = Inwii 

II dQ* , 
f>i =^V-^ V<fii- 

OWi 

Thus, dQ* /dwi acts like a force for the tpi. I now take this anal- 
ogy even further and replace it with the gradient of Q* w.r.t. the 
dependent variable ipc 



rtipi 



(where I have substituted A for trj), corresponding to 



(13) 



(14) 



For ?7 = 0, equation (13) is equivalent to the familiar equation of 
motion under the influence of the 'potential' — AQ*. Thus, if the 
time-dependence of Q* were solely due to temporal changes in the 
weights, then the energy-like quantity 



(15) 



is conserved. The frictional term proportional to ri in equation (13) 
in fact ensures that £ is not conserved but decreases, thus ultimately 
leading to the maximum of Q*, as desired. Thus, unlike the situa- 
tion for traditionally M2M, where the time-averaging of model mo- 
ments may interfere with the convergence (see the discussion in the 
last paragraph of §2), the time-averaging in my approach provides 
the damping term required for convergence. 
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3.4 Tailoring iV-body initial conditions 

As already mentioned in thie introduction, tlie M2M techinique of- 
fers a natural and powerful way to generate A''-body initial condi- 
tions with prescribed properties. There is, however, a fundamental 
difference compared to employing M2M for fitting data: the target 
values Yj now represent these prescribed properties and, unlike ob- 
served data, have no natural uncertainties. A common practice with 
model fits without known uncertainties is to simply set Uj = 1 (al- 
ternatively, setting aj = Yj alters C to measure the relative error 
and yields ST96's original method). However, this is rather unsat- 
isfactory here, as the value obtained for C then no longer provides 
insight about the goodness of fit. 

Moreover, unlike most parametric model fits, a AT-body 
model, being a Monte-Carlo representation, does have its own nat- 
ural uncertainties. This suggests that the aj should be set to the un- 
certainties expected from shot noise in the Af-body model itself. If 
this is done, C retains its interpretability: a good fit corresponds to 
C equalling to the number of constraints. While this sounds natural 
and straightforward, it introduces some subtle problems. One prob- 
lem with traditional M2M is that the time-averaging of the model 
moments intentionally reduces the shot noise, which invalidates the 
interpretability of C. 

With the alternative time-averaging of the merit function itself 
(see §3.3), this is no longer the case, but the shot noise in the N- 
body model causes temporal variations of Q* additional to those 
induced by changing the weights. This means that equation (13) 
corresponds to following a frictional trajectory (in A/'-dimensional 
y>-space) in a temporally fluctuating potential. The fluctuations are 
of the same order as the optimal value for C and prevent the algo- 
rithm to converge in the sense that oc [/, — > 0. However, the 
examples in §4 suggest that this is not a serious problem. 



3.5 Re-sampling 

The adjustment of the weights in the M2M technique may lead 
to a wide range of weights (or for Wi 7^ A^^^ to a wide range 
of Wi/wi). This potentially reduces the effective resolution of the 
iV-body model substantially and is particularly undesirable if the 
(Tj represent the uncertainties expected for a A'^-body system with 
weights following the priors. A wide range in Wi / Wi (correspond- 
ing to unequal masses for a flat prior) is also undesirable with TV- 
body initial conditions. Therefore, it is useful to re-sample the A''- 
body model from time to time during and after the adjustment pro- 
cess. This is easily done by drawing phase-space points for the new 
model from the original set {x,v)i with probability proportional 
to the relative normalised weight 

7i = w*oidM, (16) 

and subsequently setting the weights to Wi = Wi if the total weight 
is constraint to unity and Wi = WiN~^ w/c,oid otherwise. 

In this process, some trajectories of the original set will not 
be re-sampled, others get copied exactly once, yet others several 
times. In this latter case, I make the first copy a straight clone of 
the original phase-space point, but for any additional copies, I first 
randomise position and velocity as far as the underlying synmietry 
allows (for spherical symmetry, for instance, rotate them by a ran- 
dom angle about a random axis), and secondly add a small random 
velocity component. This added velocity component prevents mul- 
tiple trajectories to be on identical orbits, and aUows the model to 
explore phase-space regions of high weights. 



In order for the M2M method to still maximise the same merit 
function, one has to alter the definition of the pseudo-entropy to 

S=-^wnog'^ (17) 

with Ff the product of the factors 7i from each re-sampling so far. 

In this way, the contributions to S from each trajectory are on av- 
erage the same before and after re-sampling. However, the actual 
value for S may increase (in particular if a trajectory with 7; <C 1 
happens to be re-sampled). 



3.6 Technicalities 

The description of my M2M method, is completed by giving 
some technical details. The M2M adjustment step is taken to be 
&T = 2~^, which appears to be sufficiently short. Between these, 
the trajectories are integrated using individual adaptive time steps 
(which are required despite the fact that every trajectory is inte- 
grated on its own orbital time). While this could be implemented 
with any type of method, 1 use the traditional A/^-body block-step 
scheme with a kick-drift-kick leap-frog and a time step St = 
2Tvf^x^/\x ■ V$|/Ti with / = 1/400. In this way, trajectories 
arc automatically synchronised at M2M adjustment steps and the 
simultaneous computation of gravitational forces for many posi- 
tions allows some optimisation. 

The M2M equations are also integrated using a kick-drift-kick 
leap-frog — note that equation (14) can be integrated exactly at fixed 
Wi. In practice, e is grown slowly over 5r = 1 to its final value, but 
also Umited to e < 77 mcux{ | ?/, | } at any time. 

The AT-body model is re-sampled whenever the ratio between 
maximum and minimum w* /wi exceeds a certain threshold (4 in 
the runs of §4) and a minimum interval has elapsed since the last re- 
sampling {5t — 10 in §4). The phase-space coordinates for the fcth 
re-sampled trajectory are set to those of the ith original trajectory 
where 

Ci < 7(fe - i)< C.+i, i,k£[l,N] (18) 

with 7 = N^^ J2k "T*: mean relative normaUsed weight and 
Ci = J2k<i'yk the cumulative relative normahsed weight of the 
original model. Trajectories with 7, < 7 generate at most one copy, 
while those with 7, > 7 get copied once or more. The random com- 
ponent added to the velocities of extra copies is drawn from a nor- 
mal distribution with standard deviation 0.05 cxp(— r/10) times 
the local escape velocity (but avoiding generation of unbound tra- 
jectories). In the case of a flat prior Wi = A^ ^ ^ , which I used in §4, 
these relations simplify somewhat (in particular 7 = N~^). 

The M2M method is ideally suited for distributed-memory 
parallelisation, since the gravitational potential is fixed (no interac- 
tions between trajectories) and the evaluation of the merit function 
and its derivatives require only minimal communications. I imple- 
mented my method using the message-passing interface (MPI) and 
found the resulting code to be super-scaling: doubling the number 
of processors at fixed problem size reduces the execution time to 
less then half. This is presumably a result of the increase in to- 
tal cache, reducing the total sum of computation times, which out- 
weighs the small communication overhead. 
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Figure 1. Triaxial model: time evolution of the constraint function Cp, 
pseudo-entropy S, and rms value of the velocities of change for the M2M 
adjustment of N = 10^ particles with n = 100, e = 0.5, ?? = 0.5. 



4 APPLICATION: TAILORED JV-BODY INITIAL 
CONDITIONS 

Almost aU published iV-body simulations featuring individual 
galaxies use initially spherical dark-matter haloes with isotropic ve- 
locity distributions. This is because for these settings distribution- 
function models, on which AT-body initial conditions are usually 
based, are relatively simple to obtain. However, there is no phys- 
ical justification for these simplifications and triaxial dark haloes 
with anisotropic velocity distributions are certainly more realistic. 
Here, I apply my novel M2M method to tailor such A/^-body initial 
conditions. 



4.1 A triaxial halo model 

Let us first consider the problem of designing a triaxial equilibrium 
with prescribed shape and density profile, but without constraining 
its velocity structure. The aim is to construct a triaxial truncated 
Dehnen & McLaughlin (2005) model, which has density 



p oc 



-7/9 



4/9 



+ 1 



sech- 



n 



(19) 



with scale radius rs, truncation radius rt, and 'elliptical radius' 



2 2 2 

2 _ X y z 



(20) 



with abc = 1. For this model, the radius at which — d In p/d In g = 
2, often referred to as the scale radius for dark-matter haloes, equals 
r2 = (ll/13)^/*rs « 0.687rs. I choose rt = Wrs and axis ratios 
c/a = 0.5 and b/a = 0.7. 

A convenient way to constrain the full three-dimensional 
density distribution of the model is by means of an expansion 
in bi-orthonormal potential-density basis functions ijjnim{x) and 




100 



Figure 2. Triaxial model: density (top) and axis ratios (middle) plotted vs. 
elUptical radius for the target model (thin black lines) and the Af-body 
model just after M2M adjustment (t = 0) and after running them stan- 
dalone (under self-gravity) for 200 time units. Also plotted is the velocity 
anisotropy parameter /3 = 1 — (a'g + a'^)/2a^ (bottom). Note that the 
density at very small and large values is overestimated, an artefact of the 
density estimation procedure. 



Pnim(x)- These satisfy the Poisson equation as well as the bi- 
orthonormality and completeness conditions 



-V ^p„{x) = 4ttp„{x), 
^ Vn(a:) pn(a;') = 6{x-x'), 



(21) 
(22) 

(23) 



where n = (n, I, m). In this study, 1 use Zhao's (1996) basis set, 
whose lowest-order functions satisfy 

1 1 



V'o oc 



(|x|Va+sVa)»' 



po oc 



|a,|2-l/a(j3,jl/a_^gl/a)2+a 



(24) 



with scale radius s, and a free parameter o, which controls the den- 
sity profile. The expansion coefficients 



A„ = y^^Wjipnixj). 



(25) 



are moments of the model and of the form (3), such that the result- 
ing constraint function 

2 



A„-B„ 

<7n 



is of the form (2). Note that the calculation of the derivative 



^ 

dwi 



2^ ^"^2^" V'n(Xi) 



(26) 



(27) 



is equivalent to computing the gravitational potential due to the co- 
efficients 2{A„—B„) jda at the position Xj. This has the benefit that 
the fimctionaUty of an existing basis-function based iV-body force 
solver (dubbed 'self-consistent field code' by Hemquist & Ostriker 
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Figure 3. Similar to Fig. 1, but for the M2M adjustment of the spherical 
model with anisotropic velocities. 

1992) can be readily utilised. I use s = rg and a = 9/4 for the 
parameters of the expansion and include terms up to Umax = 20 
and Imax = 12. Since the model is forced to have triaxial symme- 
try already by the assumed underlying gravitational potential (see 
below), one does not need to constrain this symmetry. This con- 
siderably reduces the nimiber of terms in equation (26), since co- 
efficients with odd I or m can be ignored (they vanish for triaxial 
symmetry) as well as those with m, < (since An i m = i -m), 
which leaves just 588 independent constraints for A'^ = 10^ parti- 
cles. 

In order to prepare for the M2M adjustment, I sample A'^ — 
10^ positions from the target density model and evaluate the re- 
sulting A„. This is repeated many times to obtain the target val- 
ues B„ and the expected errors an from ensemble averaging. Next, 
also the velocities are sampled from the equivalent spherical (and 
velocity-isotropic) distribution-function model and scaled in each 
dimension such that the tensor virial theorem is satisfied (whereby 
correcting for unbound particles). Lastly, to achieve phase-mixing 
the resulting trajectories are integrated for several orbital times in 
the potential of the target model computed from the expansion 

$(a;) = -GMtot J2 ^" '/■'"(^)- (28) 

n 

I hoped this procedure already results in a model close to the target, 
but this was not the case at all: the resulting value for the constraint 
function (26) is - 5 x 10'' > 588. 

Finally, 1 run the M2M scheme of equations (5) and (14) with 
e = 0.5, T] = 0.5, and various values for /x, whereby integrating 
the trajectories in the target potential (28). Fig. 1 shows the time 
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Figure 4. Spherical anisotropic model; radial profiles of density and veloc- 
ity anisotropy parameter for the target model (thin black hues) and those 
measured for the AT-body model just after M2M adjustment (t = 0) and 
after rurming them standalone (under self-gravity) for 200 time units. 

evolution of Cp, S, and the rms value of Ui for an adjustment run 
with n = 100. After a quick reduction of the error (as measured by 
Cp), convergence becomes somewhat slower. Cp fluctuates with 
amplitude similar to its good-fit value of 588, as the discussion in 
§3.4 suggested, and does not reach this value. 

In order to independently assess whether the adjustment suc- 
cessfully produced a stable AT-body model with the desired proper- 
ties, I run it for 200 time units^ (corresponding to ~ 6 dynamical 
times at the scale radius) under self-gravity, whereby monitoring 
shape and density profile. This latter is done by first estimating the 
density at each particle using a kernel-estimator from its 32 near- 
est neighbours and then binning particles in density (with 5000 per 
bin) to estimate the rms radius and axis ratios (from the eigenval- 
ues of the moment-of-inertia tensor) of density shells. As evident 
from Fig. 2, the A'^-body model matches the target very well, except 
for too little flattening in the outermost regions. It appears thus, that 
the failure of convergence of Cp close to its good-fit value is caused 
by the model being too round at very large radii (well outside the 
virial radius of a CDM halo). The model appears stable: there is no 
significant change over 200 time units. 



4.2 A spherical halo model with anisotropic velocities 

Next 1 consider the problem of generating a spherical halo model 
with specified velocity anisotropy. I use the spherical version of 
the model (19) and aim to constrain Binney's anisotropy parameter 
(3=1 — {(Tg + CT^) /2(j^ to have radial profile 

(r/r l"/^ 

l3^oUr)=l3^ ^l[;/^^^,,, . (29) 

This corresponds to isotropy (/? — 0) in the very centre, and 
a slowly increasing radial anisotropy (for Poo > 0), reaching 
/3 — > /3oo at r — > oo, which describes simulated dark-matter haloes 
remarkably well (Dehnen & McLaughlin 2005). If one wants to re- 
tain the traditional M2M approach of constraining model moments, 



^ I use a unit system with = 1, G = 1, and M = 1, the total halo mass. 
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Figure 5. Similar to Figs. 1 and 3, but for M2M adjustment of triaxial model 
with anisotropic velocities. 



one must constrain the moments pcr^, pag, and pcr^ to values ob- 
tained from solving the Jeans equation. However, this latter step 
requires spherical symmetry and hence cannot be generalised to 
non-spherical systems. Instead, I directly constrain the anisotropy 
via 



Pj — Prnodclirj) 



(30) 



where j indexes radial bins with rms radius rj and measured 
anisotropy 



Hi 



(31) 



where the sums are over all particles in the jth radial bin. Since /3j 
is not a moment of the model, but a function (ratio) of moments, it 
is not of the form (3). However, the derivatives d(ij /dwu needed 
for dCp/dwi, can still be easily computed. The uncertainties (Tj 
could be estimated from the A'^-body model, but since such an es- 
timate depends on the Wi and hence adds to dC/j/dwi, this would 
complicate matters unnecessarily. Instead, I simply assume (with 
rij the number of particles in the bin) 



(1-/3) 




(32) 



100 



Figure 6. Triaxial model with anisotropic velocities: radial profiles of den- 
sity, axis ratios, and velocity anisotropy parameter for the target model (thin 
black Unes) and those measured for the AT-body model just after M2M ad- 
justment (t = 0) and after running them standalone for 200 time units. 



which is the standard deviation expected for a multivariate normal 
velocity distribution® with anisotropy j3. 

Fig. 3 shows the time evolution of the various quantities for 
an M2M run with N = 10^, rimax = 20 in Cp (only terms with 
/ = m = are considered) and 100 radial bins for Cy with f3oo = 
0.75, i.e. increasing radial anisotropy. The values for the constraint 
functions converged to their best-fit values relatively quickly and 
after r ~ 100 hardly any improvement is made. This is different 
from the situation for the triaxial halo model in the previous sub- 
section, which took much longer to reduce Cp and required much 
smaller p,. The reason for this difference is not clear, but possibly it 
is because the solution space for spherical models with anisotropic 
velocities is larger than that for triaxial models with the assumed 
axis ratios. 

The resulting radial profiles for p and (3 for the iV-body model 
provide an excellent match to the target values both just after the 
M2M adjustment and after running the model in isolation (under 
self-gravity) for 200 time units, as demonstrated in Fig. 4. 



4.3 A triaxial halo model with anisotropic velocities 

Finally, I want to generate a triaxial halo with the same triaxial den- 
sity distribution as in §4.1 but with the velocity anisotropy profile 
given by equation (29) with 0^ = 0.75, though with r replaced 
with the elliptical radius q defined in equation (20). To this end, 
I take the iV-body models generated in §4.1 as starting point for 
the M2M adjustment. Just as in the previous sub-section, the con- 
straint function now consists of two terms, constraining the density 



^ This standard deviation also depends on the ratio fl/a-^. The minimum 
occurs for = and corresponds to equation (32), while the maximum 
(arising at vanishing ag or u^) is only a factor ^4/3 larger. 
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and velocity anisotropy, respectively. Figs. 5 and 6 show, respec- 
tively, the M2M adjustment and the comparison of the final model 
with the target. Evidently, the model matches the target excellently, 
even the shape in the outermost parts, which was too round in §4.1. 



5 DISCUSSION 

The basic idea of the made-to-measure (M2M) method is to adjust 
the A'^-body weights until the model satisfies some constraints, ex- 
pressed as maximisation of a merit function. This is achieved by 
changing the weights slowly in the uphill direction of the merit 
function. While any M2M algorithm must follow this basic recipe, 
there is significant freedom in the details of how this is done. The 
purpose of this study was to improve these details compared to 
the original method as proposed by Syer & Tremaine (1996) and 
sUghtly further developed by De Lorenzi et al. (2007). 

A significant problem of this original method originates from 
the fact that the natural time scale for the adjustment (and moment- 
averaging) is some small multiple of the orbital time. Since the lat- 
ter varies substantially between orbits and, in particular, has no fi- 
nite upper limit, any finite value for the adjustment time results in 
too slow or too fast an adjustment for most orbits. I solved this 
problem by introducing a dimensionless time variable, which ef- 
fects to integrating each trajectory for the same number of orbital 
times. This is similar to Schwarzschild's method, where usually 
each orbit considered is integrated for a fixed amount of orbital 
times. In fact, the number of orbital time scales in Schwarzschild's 
method is of the same order (~ 100, depending on details, such as 
orbital symmetries) as with my M2M technique, indicating that for 
typical orbits this number is required to gather sufficient informa- 
tion. 

Note that the iterative method of Rodionov et al. (2009), men- 
tioned in the introduction, suffers from the same basic problem: 
evolving the model over some (short) time scale will bring the inner 
parts, where the dynamical time is short, much closer to equilibrim 
than the outer parts. In order to gather sufficient information about 
the dynamics in the outer parts, one would need to integrate orbits 
(or evolve the model) over a considerably longer time than is often 
practical. 

Another issue with the original M2M method is that the aver- 
aging of the model moments, required to suppress iV-body shot 
noise, interferes with the adjustment process, though apparently 
this did not lead to practical problems so far. However, if the M2M 
method is used to tailor A^-body initial conditions, the uncertainties 
entering the constraints are not observational errors but those due to 
shot noise in the iV-body model itself. In this case, time-averaging 
the model moments reduces this shot noise and renders the interpre- 
tation of the x'^-like constraint meaningless. I have overcome both 
these problems by introducing a novel adjustment algorithm which 
effects to time-averaging the merit function instead of the model 
moments and corresponds to following an orbit in AT-dimensional 
weight space with the merit function representing the potential. A 
damping term, which emerges from the averaging, guarantees that 
the maximum of the merit function will be reached. 

I also propose to (optionally) modify the merit function such 
that it automatically meets the total-mass conservation constraint. 
Finally, I propose to re-sample the A'^-body model from time to time 
during the M2M adjustment process to ( 1 ) avoid a loss of resolution 
because of unequal weights, and (2) to allow the model to explore 
phase-space in regions of high weights. This latter is achieved by 



adding a small random velocity to extra copies of trajectories, ef- 
fecting to probe another orbit close to a highly weighted one. 

Certainly, one can think of further improvements to the M2M 
method. One issue is an automatic adaption of the parameters e and 
T) to achieve optimal convergence (a technique to adapt /i such that 
the constraint function obtains a certain nvunerical value was al- 
ready proposed by ST96). The priors for the weights can be used to 
allow the Af-body model to have different mass resolution in differ- 
ent phase-space regions, which is a common technique (e.g. Zemp 
et al. 2008; Zhang & Magorrian 2008) for increasing the resolution 
in, say, the inner parts of models for dark-matter haloes. A sig- 
nificant speed-up may be achieved by starting with a relatively low 
number of particles and increasing A'^ (essentially like re-sampling) 
only later after the merit function is close to maximal. 

Finally, one would like to adapt not only the weights of the N- 
body model, but also the underlying mass distribution generating 
the gravitational potential, for instance when interpreting kinematic 
data in terms of the underlying (dark-) matter distribution. Unfortu- 
nately, changes in the orbits induced by changes of the gravitational 
potential are not straightforward to anticipate and hence to take into 
account in the adjustment process. In a spherical setting one may, 
for instance, re-scale the phase-space coordinates of every parti- 
cle such that the eccentricity, inclination, and mean radius is pre- 
served when the gravitational potential is changed. Unfortunately, 
however, something similar can no longer be done in the general, 
i.e. triaxial, case. Thus, it seems that this is a really hard problem 
and that one is forced to 'jump' from one mass model to the next 
whereby starting from the best-fit A^-body model of a 'nearby' po- 
tential. In this case, convergence may be fast, i.e. only a few ten 
orbital times, leading to significant speed-up. 

As the practical examples of the previous section demon- 
strated, my novel M2M algorithm is a powerful tool to construct N- 
body models with specified properties. One may use the method to 
explore possible stellar-dynamical equilibrivun solutions and their 
properties. For instance, the triaxial models of §4.1 exhibit a signif- 
icant radial velocity anisotropy (Fig. 2 bottom panel), even though 
the initial conditions fed to the M2M adjustment procedure were 
created from a velocity isotropic model and the velocity struc- 
ture was not constrained. This strongly suggests that radial velocity 
anisotropy is an inevitable property of (non-rotating) triaxial equi- 
librium models. This can be qualitatively understood from the in- 
evitable prevalence of box orbits, which are the only orbital family 
supporting a triaxial shape, but a more quantitative understanding 
would be desirable. 
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